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Abstract. 

This paper presents a new automatic algorithm to derive optical snow grain size (SGS) at 1 km 
resolution using Moderate Resolution Imaging Spectroradiometer (MODIS) measurements. 
Differently from previous approaches, snow grains are not assumed to be spherical but a fractal 
approach is used to account for their irregular shape. The retrieval is conceptually based on an 
analytical asymptotic radiative transfer model which predicts spectral bidirectional snow 
reflectance as a function of the grain size and ice absorption. The analytical form of solution 
leads to an explicit and fast retrieval algorithm. The time series analysis of derived SGS shows a 
good sensitivity to snow metamorphism, including melting and snow precipitation events. Pre- 
processing is performed by a Multi- Angle Implementation of Atmospheric Correction ( MAIAC) 
algorithm, which includes gridding MODIS data to 1 km resolution, water vapor retrieval, cloud 
masking and an atmospheric correction. MAIAC cloud mask (CM) is a new algorithm based on a 
time series of gridded MODIS measurements and an image-based rather than pixel-based 
processing. Extensive processing of MODIS TERRA data over Greenland shows a robust 
performance of CM algorithm in discrimination of clouds over bright snow and ice. As part of 



2 


the validation analysis, SGS derived from MODIS over selected sites in 2004 was compared to 
the microwave brightness temperature measurements of SSM\I radiometer, which is sensitive to 
the amount of liquid water in the snowpack. The comparison showed a good qualitative 
agreement, with both datasets detecting two main periods of snowmelt. Additionally, MODIS 
SGS was compared with predictions of the snow model CROCUS driven by measurements of 
the automatic whether stations of the Greenland Climate Network. We found that CROCUS 
grain size is on average a factor of two larger than MODIS-derived SGS. Overall, the agreement 
between CROCUS and MODIS results was satisfactory, in particular before and during the first 
melting period in mid- June. Following detailed time series analysis of SGS for four permanent 
sites, the paper presents SGS maps over the Greenland ice sheet for the March-September period 
of 2004. 

1. Introduction 

Snow grain is arguably the most fascinating parameter within a snowpack, because of the 
different geometries that ice grains can exhibit: from large cup shaped grains of depth hoar, to 
stellar flakes characterizing new snow, and big structures of grape-shaped rounded crystals of 
melted/refrozen snow. First observations of snowflakes were reported by Kepler (1611) and 
Descartes (Frank, 1974), the former pondering the question on why the snow crystals exhibit a 
six-fold symmetry and the latter being the first to provide a morphological description of snow 
crystals. Later on, in 1665 Robert Hooke (2003) published a volume containing sketches of 
objects observed with the microscope, including snow crystals. More recently, Wilson A. 
Bentley (1931) in 1931 performed microphotography of a large number of snow crystals and 
Ukichiro Nakaya (1954) in 1954 performed the first systematic study of snow crystals. The 
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current methods of studying the shape and size of ice crystals can be found in review by Domine 
et al. (2008). 

Grain size is a fundamental parameter for characterizing snowpack properties: for example, the 
thermodynamic state of a snowpack is related to grain size because of the metamorphism 
determined by gradients of temperature and vapor; specific snow area and snow diffusivity, 
which are primary factors in snow chemistry of gases, are related to the snow grain size; spectral 
snow albedo, playing a key role in the Earth’s radiation balance, is largely controlled by the grain 
size (e.g., Domine et al., 2008). 

Space-borne remote sensing of the Earth in the visible and shortwave infra-red (SWIR) regions 
has been used for estimating grain size at large spatial scales and moderate-to-high spatial 
resolution. In the past 10-15 years there has been a substantial body of research on retrieval of 
the snow grain size and absorptive impurities by soot or dust from satellite observations. A 
comprehensive overview of these works has recently been given by Stamnes et al. (2007). The 
proposed methods commonly rely on the fact that the absorption of pure ice is very low in the 
visible, but grows to moderate/strong in the NIR/SWIR (Wiscombe and Warren, 1980). Because 
ice is highly transparent in the visible, its reflectance has a weak dependence on the grain size 
and a high sensitivity to absorption by impurities. In the NIR region, the ice absorption increases, 
and snow reflectance becomes strongly dependent on the grain size as well as amount of soot. 
Finally in the SWIR, the ice absorption becomes dominant, and snow reflectance depends mainly 
on the grain diameter. 

The retrieval of snow grain size (SGS) relies on the model of snow reflectance as a function of 
snow properties, viewing geometry and wavelength. Commonly, Mie theory (Mie, 1908) has 
been used for this purpose, and different studies relied on the assumption of spherical particles 
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for SGS retrievals (e.g., Nolin & Dozier, 1993; Painter et al., 2003; Stamnes et al., 2007). In the 
past decade, an analytical asymptotic radiative transfer (AART) snow model was developed 
which accounts for irregular shapes of snow grains (Kokhanovsky & Zege, 2004; Kokhanovsky 
et al., 2005; Kokhanovsky, 2006). Recently, Tedesco and Kokhanovsky (2007) applied AART 
with the Koch model of fractal particles for retrieving SGS from the near-infrared measurements 
of the Moderate Resolution Imaging Spectroradiometer (MODIS). The analysis of that paper 
also identified a high sensitivity of the retrieved grain size to the errors of atmospheric 
correction. The largest impact, in particular in operational applications, is produced by errors in 
cloud detection over snow. 

In this study, we apply a new Multi- Angle Implementation of Atmospheric Correction ( MAIAC ) 
algorithm (Lyapustin and Wang, 2007) and AART snow model for the SGS retrieval over the 
Greenland ice sheet (GIS) using data collected by the MODIS instrument flying onboard the 
TERRA satellite. An improved knowledge of the temporal evolution of grain size at large spatial 
scales over Greenland is useful for refining surface melting detection techniques from other 
satellite data sets (e.g., Tedesco, 2007a) and for understanding how processes such as increased 
melting (e.g., Tedesco, 2007b) affect the surface energy balance and vertical stratigraphy and 
precipitation events (in both timing and extent) over the GIS. Because of the paucity of grain size 
measurements over the GIS and of the dependency of these measurements on the operator’s 
ability and skills, results derived from MODIS data are compared with those obtained with a 
snow model driven with surface meteorological data recorded on ground by the Greenland 
Climate Network (GC-Net) automatic stations. 

The paper is structured as follows: in Section 2 we discuss MAIAC cloud detection and 
atmospheric correction, and provide detail of the AART snow model. The algorithm of SGS 
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retrieval is described in section 3. In Section 4 we focus on the results obtained for selected sites 
over the GIS, and on the comparison between satellite-derived SGS and those obtained with a 
snow model. A qualitative analysis between MODIS SGS and microwave data measured by the 
Special Sensor Microwave Imaging (SSM/I) radiometer is also reported. We discuss the results 
in section 5 and dedicate the last Section to the summary. 

2. Method Description 

MAIAC is a new algorithm which uses a time series of MODIS measurements and combines 
image- and pixel-level processing to detect clouds and perform simultaneous retrievals of aerosol 
properties and surface bidirectional reflectance and albedo (Lyapustin and Wang, 2007). The 
algorithm is generic although the aerosol retrievals are currently not performed over snow. In 
this work, we are using a reduced MAIAC processing, which includes cloud mask, water vapor 
retrievals and a simplified form of atmospheric correction described in sections 2. 1-2.2. 

2.1 Cloud Mask Algorithm 

The central idea of MAIAC cloud mask (CM) algorithm is to build and dynamically update a 
reference clear-skies image of the surface ( refcm ), which is further used as a comparison target 
for cloud detection. With a high frequency of global MODIS observations, the land surface can 
be considered as a static or slowly changing background contrary to ephemeral clouds. In cloud- 
free conditions, the consecutive images of the same surface area have a well-reproducible spatial 
pattern and a high spatial correlation. Clouds introduce distortions in this pattern usually in a 
random way reducing the correlation, which can be detected by covariance analysis. The 
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covariance is a metric showing how well the two images X and Y correlate over an area of NxN 
pixels, 
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Because covariance removes the average component of the signals, this metric is equally 
successful over the dark and bright surfaces and in both clear and hazy conditions if the surface 
spatial variability is still detectable from space. 


The MAIAC algorithm starts with accumulation of 6-16 days of MODIS data in the processing 
Queue. The size of the Queue depends on the frequency of MODIS observations over a given 
location, which is defined by the latitude. MAIAC uses 16 days for the Earth equator where 
MODIS makes one observation per day, and 6 days for the poles, which feature several orbits per 
day. For the time series analysis, MODIS data are gridded to 1 km resolution, creating Tiles of 
data. 


The covariance analysis is applied for the areas of 25x25 km , called blocks, using band 1 (0.645 
pm) data. When the cloud-free conditions for a given block are found, the image is copied to the 
reference image ( refcm ), which is further used to compute covariance with the new MODIS 
measurements. Each time the new measurements are identified as cloud-free for a given block, 
the refcm is updated thus dynamically adapting to seasonal variations of the landcover. In order 
to account for the effects related to scan angle variation, e.g. pixel size growth, surface BRF 
effect or reduction of contrast at higher view zenith angles (VZA), two reference clear-skies 
images are maintained by the algorithm, refcml for the view zenith angle range VZA«0-45° and 
refcm2 for VZA=45°-60°. In addition to refcm, the maximal value (rl max ) and the variance (cp) of 
reflectance in band 1 as well as the brightness temperature contrast (ABT=BT / „ av -BT„„„) for each 
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surface block are also stored in the memory. Analysis of MODIS data shows that thermal 
contrast (ABT) is a rather stable metric of a given land area in clear conditions. In partially 
cloudy conditions, the contrast increases because BT mi „ is usually lower over clouds. 

The new CM algorithm has an internal surface classifier, producing a dynamic land-water-snow 
( LWS) mask, and a surface change mask. These are an integral part of MAI AC guiding the cloud 
masking when the surface changes rapidly as a result of fires, floods or snow fall/ablation. 

The final values of the MAIAC cloud mask are clear (CM CLEAR, CMCLEARWATER, 
CM CLEAR SNOW), indicating surface type as well, possibly cloudy (CM PCLOUD), and 
confidently cloudy (CM CLOUD). The covariance component of MAIAC algorithm, which 
offers a direct way to identify clear conditions, renders another commonly used value of cloud 
mask - “possibly clear” - redundant. The full description of the algorithm, including the logic of 
cloud detection at a pixel level, is given in (Lyapustin et al., 2008). 

The current version of CM algorithm was developed for generic land applications. It may need 
further enhancements for Arctic applications, such as discrimination of the sea ice. Nevertheless, 
we found that it has a good performance over the snow detecting significantly more clear snow 
pixels than the operational MODIS cloud mask algorithm MOD35 (Lyapustin et ah, 2008). 

An example of MAIAC cloud mask for Thule located on the north-west of Greenland (77°30’0N, 
69°10’0W), is shown in Figure 1. The image shows 6 consecutive daytime MODIS TERRA 
gridded images for two days 1 12-113, 2004, for the area of 150x150 km . The two left columns 
show the TERRA RGB images. They are differently normalized to help distinguish clouds from 
the bright snow. The central feature in the image shows a fiord covered by snow/ice, which melts 


in the summer. 
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The right column shows the derived cloud mask, which uses the following legend: red - cloud; 
white, blue and light blue - clear over snow, land and water, respectively. A reproducible spatial 
pattern helps to identify cloud-free blocks of the surface. The top three images represent the 
clear-skies conditions. The second image is obtained near the edge of scan and has an apparently 
coarser resolution than images 1 and 3, acquired close to nadir. Clouds are seen in the images 4-6 
as distortions to the spatial pattern of reflectance. For the most part, they are captured by the 
cloud mask algorithm. 

2.2 Atmospheric Correction of MODIS Data 

The aerosol concentration over Greenland is generally low (Tomasi et al., 2007), except events 
of transported biomass burning smoke or pollutions from Northern America or Eurasia. Because 
of high surface reflectivity and cloud contamination, reliable remote sensing methods for aerosol 
retrievals over this region of the world do not exist. For these reasons, we only consider effects 
of Rayleigh scattering and molecular absorption, including that of water vapor and ozone. 

Table 1 shows the effective center wavelengths of MODIS TERRA 500 m bands, the in-band 
optical thickness of absorption by well-mixed gases, and of water vapor absorption at the column 
amount of 1 cm. The effective band center wavelength is defined as a wavelength for which 
monochromatic and narrow-band direct vertical transmittances of the aerosol-free atmosphere 
are equal: 

exp {-t r (X c )} = \f x exp{-r*(A )}h x dA / J F x h x dA , (2) 

A , 1 / AA 

where F x is extraterrestrial spectral solar irradiance, h x is relative spectral response (RSR) 
function of MODIS TERRA bands, and is Rayleigh optical thickness computed with model of 
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Bodhaine et al. (1999). The optical thickness of the in-band gaseous or water vapor absorption is 
defined similarly. The optical thickness of absorbing gases is calculated for the carbon dioxide 
concentration of 380 ppm, and concentration of other gases corresponding to the Arctic 
atmospheric model (Kneizys et al., 1996). The calculations included absorption of 6 major 
atmospheric gases (H2O, CO2, CH4, NO2, CO, N2O) with line parameters from HITRAN-2000 
(Rothman et al., 2003) database using the Voigt vertical profile, and the Atmospheric 
Environmental Research continuum absorption model (Mlawer et al., 2006). 

The top-of-atmosphere (TOA) reflectance is modeled in this work using a modified Lambertian 
formulation, and assuming that most of the ozone absorption takes place above the molecular 
atmosphere: 


R aa Oo =( R 'L ( Ao >M><P) + 




)T°J n \p„p) 


( 3 ) 


where R A is atmospheric path reflectance, T is the total upward and downward atmospheric 
transmittance, 5 is spherical albedo of the atmosphere, p and q are the surface bidirectional 
reflectance factor (BRF) and albedo, respectively. The subscript indicates that the functions are 
weighted with respective RSR of MODIS channels: 


R AA =\ F A R A h i d ^l\ F A h x dA . (4) 

Because atmospheric optical thickness over Greenland is low, correction using Lambertian 
equation (3) is rather accurate even at relatively low sun zenith angle (Lyapustin, 1999). The 
atmospheric correction is based on the look-up table computed using radiative transfer code 
SHARM (Lyapustin, 2005) and the Interpolation and Profde Correction (IPC) method 
(Lyapustin, 2003) for narrow-band and shortwave broadband calculations. 
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After generating cloud mask, MAI AC retrieves column water vapor with a relative accuracy of 5- 
10% from MODIS NIR channels 17-19 located in the water vapor absorption band 0.94 pm 
using a band ratio technique (Lyapustin and Wang, 2007). The retrieved values over Greenland 
for the period covering spring to fall of 2004 agree well with known climatology showing the 
average values of 0.2-0.4 cm for the spring and autumn, and 0.4-0. 8 cm for the summer period 
(Kneizys et al., 1996). 

Because of significant surface height variation over Greenland, the height/pressure correction is 
implemented for functions R A , T and s in the visible and near infrared bands 1-4 as described in 
(Lyapustin and Wang, 2007). 

2.3 Snow Reflectance Model 

In this work, we are using the AART model of snow reflectance developed by Kokhanovsky and 
Zege (2004). This model predicts snow reflectance using an asymptotic analytical solution of the 
radiative transfer problem in semi-infinite media with low absorption (the so-called Milne 
problem), which is a good approximation for snow in the visible spectral range: 

P = R o (Ao > A, <P) exp(-^(// 0 » A, <P \ ) • (5) 

Here, y = 4/zj / A , % is imaginary part of refractive index of ice, X is the wavelength, and 
d=6(y)l(S) is the effective grain size defined by the ratio of the average volume to the average 

surface area of grains. Ro is a solution of radiative transfer equation for snow with zero 
absorption. For simplicity, we assume that Ro does not depend on snow grain size because the 
wavelengths at consideration (0.4—2. 1 pm) are much smaller than parameter d, which is known 


11 


to vary from about 50 jam for fresh snow to 1-1.5 mm in pre-melting conditions (Wiscombe and 
Warren, 1980) to several millimeters during snowmelt. The function A relates to the photon’s 
escape probability from the media. For the fractal ice particles, it can be approximated as 
follows: 

AMo = 0.66(1 + 2ju 0 )(1 + 2 ju) / R 0 (// 0 , /u, (p ) . (6) 

Comparison of the described analytical model with the accurate radiative transfer solution 
showed the accuracy of better than about 3% for zenith angles below «78°. A local validation of 
the model in the spectral range of 0.535-2.21 pm with field measurements has also demonstrated 
a good angular and spectral agreement (Kokhanovsky et al, 2005). 

The described snow reflectance model has several limitations. First, as we mentioned, the zenith 
angles cannot be too high, which may limit retrievals over polar regions on winter, early spring 
and late fall when solar zenith angle is low. Second, comparisons with measurements showed 
that the model is not applicable near the forward scattering (glint) direction (tp~0°). Third, the 

model was derived with the assumption of low absorption, J3 = *Jyd « 1. This assumption holds 
well for MODIS vis-NIR bands 1-5, but in the shortwave infrared region for bands 6 (1.6 pm) 
and 7 (2.1 pm) parameter (3 may be as high as 0.5-1 during snowmelt. On the other hand, 
ground-based measurements of the bidirectional reflectance of snow (Kokhanovsky et al, 2005) 
and processing of AATSR data (Kokhanovsky and Schreirer, 2007) showed that the model can 
be used in the SWIR region as well. 

In general, the snow model should include the absorption of light by liquid water in snow (Green 
et al., 2002), and also by absorbing impurities, such as soot from biomass burning and urban- 
industrial aerosols (Warren, 2007), or transported dust. In this work, however, we did not find 
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the evidence from MODIS data of measurable soot contamination over the Greenland in 2004, at 
least above sensitivity of our model. Our analysis also showed that adding water absorption has 
almost no effect in the MODIS bands used for the grain size retrievals because the imaginary 
part of the refractive index is nearly the same. For these reasons, we decided not to consider 
these sources of absorption further. 

An obvious advantage offered by this model is an explicit dependence on the grain size and snow 
refractive index. This model offers a straightforward and very effective inversion scheme, which 
would only require tabulating function Ro as a function of angles. For this work, the LUT of 
function Ro was computed with the radiative transfer code SHARM (Lyapustin, 2005) using the 
scattering function for random fractal crystals (Mischenko et al., 1999). Mischenko provides 
2000 Legendre expansion coefficients. We found that SHARM solution achieved convergence 
with an accuracy of about 0.1 -0.2% when we used 1024 coefficients. 

The last row of Table 1 shows the band-average imaginary part of an effective refractive index of 
ice xT > based on the latest compilation of Warren and Brandt (2008), which is computed from 
equation 

exp {~A{ju 0 , ju, <P)- ] J 4 ^ d ) = J exp(-d(/v 0 , //, < p) J^5^ )F 2 h /i dA j^F^h^dA. (7) 

The effective value accounts for spectral variations of x" e (^) within the band-pass interval, for 
variations of function A=0.4-2.2 with the MODIS viewing geometry, and for changes of diameter 
during snowmelt. We found that using effective values improves spectral agreement between the 
model and the MODIS measurements in the NIR and SWIR bands. 


3. Method of Grain Size Retrievals 
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The developed algorithm has several processing stages: a) gridding MODIS level IB geolocated 
and calibrated data to 1 km resolution; b ) retrieving water vapor; c) cloud masking; d) 
performing atmospheric correction; e ) retrieving snow properties. 

The snow grain size can be derived from satellite data using either a single channel (e.g. Nolin 
and Dozier, 1993; Stamnes et ah, 2007), or often used band-ratio algorithms (e.g., Scambos et 
al., 2007). The first approach requires an accurate knowledge of the bidirectional reflectance Ro. 
We analyzed the accuracy of the modeled function Ro experimentally by comparing it with the 
atmospherically corrected MODIS TERRA reflectance for a full range of viewing geometries 
and conditions during 7 months (March - September) of 2004 over Greenland. Figure 2 shows a 
comparison for the Blue (B3) band for the cloud-free observations, as detected by MAIAC cloud 
mask for the Saddle site (66.0006N, 44.5014W). Figure 2a shows the results for all view angles 
with geometry plotted in the bottom Figure 2c. The spikes of the function Ro correspond to the 
glint directions where the AART model is not valid. It is interesting that MODIS measurements 
do not show significant increase of reflectance in the forward scattering directions. This result is 
typical not only for this particular site, but for the full area of Greenland. Figure 2b shows a 
reduced set of measurements where azimuthal angles (p<40° were filtered. Although the 
agreement in this case is significantly better, the MODIS reflectance still shows less anisotropy 
than theoretical predictions. The main reason for this disagreement is a macroscopic surface 
roughness caused by the local terrain variations (slope and height variations), sastrugis, etc. 
There is a significant amount of such variability inside 1 km grid cell which flattens the BRF 
pattern of reflectance, and which is not taken into account by the plane-parallel radiative transfer 
model used to compute Ro (e.g., Hudson et al., 2006; Hudson and Warren, 2007). 
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For this reason, we selected the band-ratio algorithm where the role of function Ro is reduced to 


the second-order effect manifested through function A. Using two channels with different ice 
absorption, it is easy to obtain from Eq. (5): 


Because of Ro, there is some uncertainty in the magnitude of function A in equation ( 8 ), which 
may introduce some bias in the value of SGS. This function, however, depends only on the 
viewing geometry, so the possible bias should not affect the spatial variations of derived SGS. 
The calibration of the retrieved snow grain size with ground-based measurements can be used to 
improve the model of function A. Some preliminary results on such experimental calibration 
have been reported by Kokhanovsky (2006). 

The example of retrieved snow grain size using the ratio R1.24/R0.65 is shown in the last column of 
Figure 1. Gridding allows an easy comparison between successive observations of the same area. 
The spatial pattern of the snow grain size is reproduced consistently for different observations on 
a single day and on successive days. The residual dependence on the VZA and on azimuth is 
small. In the retrievals, the forward scattering directions (p<40° were excluded. 

A time series of retrieved grain size for the Saddle site is shown in Figure 3. Here, we used 
several different ratios created from bands 1 (0.645 pm), 2 (0.87 pm), and 5 (1.24 pm). The 
ratios R 5 /R 1 and R 5 /R 2 produce essentially the same results, showing very stable retrievals when 
no melting occurs, and a similar growth of SGS during the two periods of snowmelt. Another 
ratio R 2 /R 1 (red line) used by Scambos et al., (2007) has a good sensitivity during snowmelt 



( 8 ) 


periods, but produces unstable results for the other periods with values often below 50 pm. Our 
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results suggest that the difference in the ice absorption between the Red and NIR bands may be 
too low for the reliable SGS retrievals when the grain size is small and no melting takes place. 

The blue line shows SGS retrieved using Eq. (5) from a single band 6 (1.64 pm) where the ice 
absorption is high. A similar approach was used for SGS retrievals from the Global Imager (GLI) 
(Stamnes et al., 2007). In complete agreement with results from the above reference, this 
retrieval shows smaller values of SGS as compared to shorter wavelengths. The higher values of 
the retrieved grain size at shorter wavelengths are usually explained by a larger depth of 
penetration of light into the icepack (e.g. Li et al., 2001) and growth of the effective grain size 
with depth (Grenfell et al., 1994; Aoki et al., 2000). 

Further, the ratio R5/R1 will be used as a baseline for our retrievals. In order to assess the 
sensitivity of the method we propose, the retrievals were conducted for different sites in the 
inland Greenland and compared with outputs of a snow energy balance model and the trends of 
space-borne microwave brightness temperature, which is described in the next section. 
Independently, we studied the behavior of SGS in the coastal regions of Greenland characterized 
by persisting melting and variations due to re-freeze and/or new snowfall. Figure 4 shows 
retrievals for the glacial area north-east of Thule shown by a red circle in the top-left image of 
Figure 4. There are several pronounced periods of systematic growth of the SGS, followed by a 
sharp drop. Some of these are selected by ovals in the figure. We have looked at periods 1 and 2 
in detail trying to identify whether this is an uncertainty of the solution or it relates to a change of 
the state of snow. In both cases, the drop of SGS was preceded by a period of cloudiness, and it 
correlates with decrease of the measured SSM/I brightness temperature. The spatial RGB images 
are shown next to the time series of SGS. We cannot provide any detail for the first period (days 
134-136) except clearing of the water from the sea ice. However, there was an evident fresh 
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snowfall during the extended 6-day period of cloudiness for the second period (days 211-217). 
On day 21 1 we see snow-free land along the melted water of the fiord. After 6 cloudy days, most 
of the bare land is covered by a new snow on day 217, with visible SGS decrease. This example 
shows that the developed method is sufficiently sensitive to detect snow metamorphism, 
associated with melting, re-freeze and fresh snowfall. 

4. Comparisons with snow model outputs and microwave data 

In this section, grain size values retrieved with MODIS are evaluated vs. those simulated with a 
snow physical model, CROCUS (Brun et al., 1989 and 1992). The need of comparing satellite- 
derived grain size with those obtained from a physical model is dictated by the lack of ground 
measurements, especially over the Greenland ice sheet. Even when they exist, these 
measurements are limited in space and time and therefore are limited in their usefulness for the 
evaluation of the algorithm proposed here. 

Snow grain in CROCUS is defined in terms of optical size for albedo computation and physical 
size for stability assessment. The state of snow grains is also described through sphericity and 
dendricity (Brun et ah, 1992). For the fresh snow, an initial value of 100 pm is set to the optical 
grain size while the initial values of dendricity and sphericity are obtained according to the 
surface wind. Then the snow metamorphoses and when the dendricity reaches 0 (from an initial 
value of -1), all the branched crystals have almost disappeared, and a size of 300 pm is applied. 
Grain size evolution is then computed according to the forcing data. In the following, forcing to 
the CROCUS model is derived from measurements of the automatic whether stations of the 


Greenland Climate Network (GC-Net, Steffen et ah, 1996). 
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Beside the results from CROCUS, a qualitative analysis of the MODIS-retrieved grain size is 
also performed using data of the Special Sensor Microwave Imaging (SSM/I) radiometer. The 
microwave brightness temperature at 19.35 GHz, horizontal polarization, recorded by SSM/I is 
plotted for the cells containing the test site at resolution of 25 km. Microwave brightness 
temperature at 19.35 GHz is sensitive to the liquid water content within the snowpack (e.g., 
Tedesco et al., 2007). When snow melts, snow grains tend to bound to each others (constructive 
metamorphism) with the grain size of melted and/or refrozen snow generally greater than that of 
unthawed snow. 

The comparison between MODIS-retrieved and CROCUS simulated SGS values, and brightness 
temperature was performed over four selected test sites: Dye-2 (66.48 ION, 46.2800W), Saddle 
(66.0006N, 44.5014W), South Dome (63.1489N, 44.8167W) and NASA SE (66.4797N, 
42.5002W). The results are shown in Figure 5. The MODIS results are averaged over an area of 
25x25 km 2 for an easy comparison with SSM\I measurements. Note that CROCUS SGS is 
reduced by a factor of two. In this case, there is a good general agreement between CROCUS 
values and MODIS retrievals for the off-melting season and the first period of strong thawing 
(June, days 170-176). The best agreement for this period is observed for NASA SE and Saddle 
sites. For the Dye-2 and South-Dome sites, the maximum of melting is shifted by several days. 

Overall, analysis of MODIS results suggest two well-expressed melting periods during the 
summer of 2004, with the second one occurring in August (days 220-226). One can see that 
CROCUS overestimates the length of melting period by almost 10 days over Saddle and South- 
Dome, but it fails to predict a strong melting event over site Dye-2. The SSM\I brightness 
temperature for this period of time supports the results obtained with MODIS. 
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Between June and August, CROCUS simulations predict additional intermediate high amplitude 
thawing periods for sites NASA-SE, Saddle and South-Dome, which are not detected with 
MODIS. Some of the disagreement between MODIS and CROCUS might be caused by 
cloudiness. Because of the inevitable errors in cloud detection, when cloudy pixels are masked as 
clear, the developed algorithm may produce lower SGS values in cloudy conditions. On the other 
hand, SSM\I data, which are not affected by cloudiness, do not support the predicted melting by 
CROCUS. 

5. Seasonal melting over Greenland in 2004 

Monitoring SGS over GIS can support the development of algorithms for snowmelt detection. 
The snow grain size retrieved from MODIS TERRA data for the melting period of 2004 over 
GIS is shown in Figure 6. Each image is obtained as a 3-day composite to cover gaps caused by 
clouds. Despite some inner areas remain un-initialized with our retrievals (shown by black 
color), the images provide a sufficiently good coverage and temporal resolution to detect short 
periods of melting. 

The first image shows SGS for the middle of April, 2004. Note that we don’t currently 
discriminate sea ice which has similar characteristics to snow of the inner regions of GIS. In 
May- June, melting takes place in the peripheral area of GIS and adjacent sea ice. The red color 
on the border of GIS indicates high SGS values. On the other hand, a relatively coarse resolution 
of our data (1 km) may result in undetected presence of pools of pure water cleared from sea ice 
along the shores, which will also be interpreted as higher SGS values. 

Figures for the last week of June (days 169-175) show strong melting in the southern inner 
region of Greenland, as well as along the northern and western perimeter of GIS. The melting 
largely abates by day 178. In July, melting continues mainly in the north of Greenland. A second 
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period of melting begins in the middle of August. The strip of low values in the middle of 
southern GIS on day 223 is probably an artifact as this region was continuously covered by 
clouds. 

We also studied the capabilities of our algorithm to produce spatially distributed SGS. A first 
step in this direction is to compare the maps of microwave brightness temperature with those of 
SGS. Figure 7 shows brightness temperature images of SSM\I for selected days of 2004. Low 
brightness temperature values are characteristics of dry snow whereas relatively high values are 
due to the presence of wet snow. From the visual comparison between microwave data and 
MODIS SGS, we observe that high brightness temperature values (melting snow) correspond to 
areas where grain size values are high. This is consistently true for all considered days. 


6. Conclusions 

Snow grain size is an important parameter, knowledge of which is required in various 
disciplines, from snow chemistry of gases to the Earth’s radiation balance. The grain size, along 
with snow temperature and soot concentration, was an operational product of the Global Imager 
aboard Advanced Earth Observing Satellite-II (Hori et al., 2007) during its mission in 2003. In 
this paper, we described a new algorithm with operational capabilities to derive snow grain size 
from MODIS data over Greenland. 

The new algorithm uses the framework of MAI AC processing to detect clouds, retrieve water 
vapor and perform atmospheric correction of MODIS measurements. Cloud discrimination over 
bright snow and ice surfaces is a very difficult problem, and robust performance of MAI AC CM 
algorithm is a key to successful retrieval of snow properties. 
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The SGS retrieval is based on the AART model, which predicts analytically the spectral and 
directional snow reflectance as an explicit function of snow grain size and absorption within the 
snowpack. We took advantage of the analytical form of solution to develop an explicit and fast 
SGS retrieval algorithm using the ratio of MODIS measurements in bands 5 (1.24 pm) and 1 
(0.645 pm). Additional retrievals in the SWIR bands 6 (1.6 pm) and 7 (2.11 pm) with higher 
show absorption and lower light penetration depth can be used to evaluate the vertical profile of 
SGS (Li et al., 2001). The consistency of solution was demonstrated over the full range of the 
MODIS TERRA viewing geometries, excluding the forward scattering (glint) range of azimuthal 
angles (cp<40°), where the AART model is not valid. We also demonstrated the consistency of 
the time series of retrieved SGS and a good sensitivity of the algorithm to detect snowmelt and at 
least some of the snow precipitation events. The results of retrievals were compared for the year 
of 2004 with the microwave brightness temperature of SSM\I radiometer which is sensitive to 
amount of liquid water in the snow. The Greenland-scale comparisons show a good agreement in 
detection of two main melting periods in 2004, in the middle of June and in August. We believe 
that MODIS-based SGS retrievals may have a better sensitivity to the changes in snow 
metamorphism than the microwave measurements, in part due to significantly better spatial 
resolution (1 km vs 25 km), though temporal resolution might be affected by clouds presence. 

MODIS retrievals were also compared with the snow model CROCUS driven with surface 
meteorological data recorded on ground by the Greenland Climate Network automatic stations 
for four different selected sites. Overall, the comparison was satisfactory, although the grain size 
of CROCUS is found to be on average a factor of two larger, and the number of predicted 
periods of thaw higher. There are certain differences in the magnitude and the longevity of the 
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main melting periods as well. The independent SSM\I results agree much better with MODIS 
retrievals than with CROCUS prediction results. 

The performed comparison with SSM\I and CROCUS data, as well as good quantitative 
agreement with validated and published GLI results (Hori et al., 2007; Aoki et al., 2007) serve as 
an indirect validation of the described method. In the near future, we plan to extend validation to 
comparison with ground-based albedo measurements as well as available in situ grain size 
measurements. 

The current algorithm is fast and can be run operationally. On the other hand, our processing 
model, although ensuring robust performance over Greenland, needs further development. We 
plan to include the light absorption by snow impurities, such as soot, in order to adapt the 
algorithm for the regions with high urban impact on snow. An enhancement is needed to 
discriminate sea ice based on MODIS data, and fine-tune MAIAC cloud mask for the Arctic 
environment. We demonstrated a sensitivity of method to detect fresh snowfall, but further work 
is needed to get to the point of automated classification of such events. One important part of 
research is developing a theoretical model to account for the surface roughness leading to 
reduced anisotropy of snow reflectance (function Ro) as compared to the currently used plane- 
parallel modeling approach. In this sense, MODIS is an excellent source of data with high 
frequency of observations over polar regions and a wide angular coverage. 
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BandN 

B1 

B2 

B3 

B4 

B5 

B6 

B7 

A c , pm 

0.6449 

0.8556 

0.4655 

0.5535 

1.2419 

1.6290 

2.1131 

T g 

1.32e-3 

1.82e-5 

2.06e-3 

4.60e-4 

2.90e-3 

1. Ole-2 

1.80e-2 

r w , 1 cm 

3.62e-3 

5.49e-3 

5.20e-5 

3.63e-4 

3.62e-3 

0.87e-3 

1.63e-2 

y ice 
A / e 

1.25e-8 

2.32e-7 

1.05e-9 

3.22e-9 

1.20e-5 

2.41e-4 

(5.3-6.8)e-4 


Table 1 . The center wavelength of MODIS TERRA land 500 m channels, the in-band 
absorption optical thickness of well-mixed gases and of 1 cm of water vapor, and effective 
refractive index of ice. 
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Figure 1 . Illustration of MAI AC cloud mask for two consecutive days 112-113 of 2004 of 
MODIS TERRA observations over Thule, Greenland (the area of 150x150 km 2 ). The left 
two columns show the RGB image of gridded MODIS LIB data. They are normalized 
differently to help distinguish clouds over the snow. The third column shows MAI AC cloud 
mask, which uses the following legend: red - cloud; white, blue and light blue - clear over 
snow, land and water, respectively. The last column shows the retrieved snow grain size (see 
section 3). 
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Figure 2. Comparison of function Ro (solid line) with atmospherically corrected MODIS 
TERRA reflectance in the Blue band (0.47 pm) for the Saddle site. The top figure (a) shows 
the full set of clear-skies measurements. The bottom figure (b) shows the reduced dataset 
with forward scattering directions (cp<40°) filtered. The bottom plot (c) shows the geometry 
of observations for the full set of measurements. The smooth solid line and dashed line with 
circles show the cosine of solar and view zenith angles, respectively. The rapidly changing 
solid line shows the relative azimuth. 
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Figure 3. A time series of the derived snow grain size over Saddle, Greenland, for the period of 
March - September 2004. The retrievals used the following channels: B1 (0.645 pm), B2 
(0.87 pm), B5 (1.24 pm), B6 (1.64 pm). The x-axis shows the number of MODIS TERRA 
consecutive observation. 
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Figure 4. (Top) A time series of the snow grain size and brightness temperature for the area 
shown by a red circle near Thule, Greenland. (Bottom) MODIS TERRA RGB images, cloud 
mask, and retrieved snow grain size for the periods 1 (days 134-136) and 2 (days 211-217). 
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Figure 5. Comparison of the time series of snow grain size and SSM\I brightness temperature for 
four sites of Greenland. The gray line show SSM\I brightness temperature, the solid line 
shows CROCUS SGS reduced by a factor of two, and circles connected by a dashed line 


show MODIS-retrieved SGS. 
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Figure 6. Snow grain size (mm) retrieved from MODIS TERRA data for the melting period of 
2004 over Greenland. Each image is obtained as a 3 -day composite to cover gaps caused by 
clouds. 
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Figure 7. SSM\I brightness temperature over Greenland for 2004. 


